Reference site numbering¶
So far all the preceding examples have involved analysis of proteins with sequential integer site numbers (e.g., 1, 2, 3, …). This sequential integer numbering can also start at sites other than one (e.g., 331, 332, 333, …).
However, often proteins are numbered in relation to alignment to sequential integer numbering of some reference homolog. For instance, SARS-CoV-2 spike sequences are usually numbered in relation to the Wuhan-Hu-1 reference. But if the protein has accumulated insertions or deletions relative to the reference, there may be missing sites (gaps) or insertions (typically numbered as 214, 214a, 214b, etc).
This notebook shows how to perform an analysis and analyze the results with non-sequential reference site numbering. The advantage of doing this is that the results can directly be visualized/analyzed using the conventional site numbering scheme.
It does this by analyzing deep mutational scanning of the SARS-CoV-2 Omicron BA.1 spike, which has several indels relative to the Wuhan-Hu-1 numbering reference.
Fit and visualize a model with non-integer site numbering¶
First, import Python modules:
[1]:
import polyclonal
import pandas as pd
Now read the data to fit:
[2]:
# read data w `na_filter=None` so empty aa_substitutions read as such rather than NA
data_to_fit = pd.read_csv(
"Lib-2_2022-06-22_thaw-1_LyCoV-1404_1_prob_escape.csv",
na_filter=None,
)
data_to_fit.head()
[2]:
| antibody_concentration | barcode | aa_substitutions_sequential | aa_substitutions_reference | n_aa_substitutions | prob_escape | prob_escape_uncensored | no-antibody_count | |
|---|---|---|---|---|---|---|---|---|
| 0 | 0.654 | AAGAGCTTACTCTCGA | S443P S1167P | S446P S1170P | 2 | 0.9004 | 0.9004 | 2890 |
| 1 | 0.654 | ATAAGATAGATTTAGG | K441T A519S | K444T A522S | 2 | 0.8006 | 0.8006 | 2217 |
| 2 | 0.654 | GATCGAGTGTGTAGCA | F152T K441E | F157T K444E | 2 | 0.5797 | 0.5797 | 2966 |
| 3 | 0.654 | CCAAACGGTATGATGA | Q14H V67P D212H N445T I1224M | Q14H V67P D215H N448T I1227M | 5 | 0.4125 | 0.4125 | 4131 |
| 4 | 0.654 | TTACTGTGCAACCCAA | F163L N445S | F168L N448S | 2 | 0.3600 | 0.3600 | 4275 |
Notice how these data have sites numbered in two different schemes:
aa_substitutions_sequential: sequential integer numbering (1, 2, 3, …) of the protein used in the experiment. You could just analyze the mutations this way, but then the results will not be in standard reference numbering scheme.aa_substitutions_reference: the referencey based site numbering, which skips sites with indels and has some non-numeric sites where there are insertions (eg,214a), as for example in the variant shown below:
[3]:
data_to_fit.query("aa_substitutions_reference.str.contains('a')").head(n=1)
[3]:
| antibody_concentration | barcode | aa_substitutions_sequential | aa_substitutions_reference | n_aa_substitutions | prob_escape | prob_escape_uncensored | no-antibody_count | |
|---|---|---|---|---|---|---|---|---|
| 32 | 0.654 | CTAGATAAATCCCTGC | E209A K441N | E214aA K444N | 2 | 0.3965 | 0.3965 | 2214 |
Here we will use the reference based amino-acid substitutions, so create a column with the name used by Polyclonal (the “aa_substitutions” column) that uses this numbering scheme:
[4]:
data_to_fit["aa_substitutions"] = data_to_fit["aa_substitutions_reference"]
Importantly, in order to use reference based numbering that is not sequential integer, you also have to provide a list of the sites in order. The reason is that otherwise it’s not possible for Polyclonal to figure out for instance if sites are just missing from the data are are actually deletions.
Here we read a data frame that maps sequential to reference site numbering, and then use that to extract the list of reference sites:
[5]:
site_numbering_map = pd.read_csv("BA.1_site_numbering_map.csv")
display(site_numbering_map.head())
print("Note how some reference sites differ from sequential ones due to indels:")
display(
site_numbering_map[
site_numbering_map["sequential_site"].astype(str)
!= site_numbering_map["reference_site"]
].head()
)
sites = site_numbering_map["reference_site"].tolist()
| sequential_site | reference_site | |
|---|---|---|
| 0 | 1 | 1 |
| 1 | 2 | 2 |
| 2 | 3 | 3 |
| 3 | 4 | 4 |
| 4 | 5 | 5 |
Note how some reference sites differ from sequential ones due to indels:
| sequential_site | reference_site | |
|---|---|---|
| 68 | 69 | 71 |
| 69 | 70 | 72 |
| 70 | 71 | 73 |
| 71 | 72 | 74 |
| 72 | 73 | 75 |
Now initialize and fit the Polyclonal model, but pass the sites argument so we can use these non-sequential-integer reference sites. (If you don’t pass the sites argument, sites are assumed to be sequential integer):
[6]:
model = polyclonal.Polyclonal(
# `polyclonal` expects the concentration column to be named "concentration"
data_to_fit=data_to_fit.rename(columns={"antibody_concentration": "concentration"}),
n_epitopes=1,
alphabet=polyclonal.AAS_WITHSTOP_WITHGAP,
sites=sites,
)
Note how the model has its sequential_integer_sites attribute set to False:
[7]:
assert set(model.sites) == set(sites)
assert model.sequential_integer_sites is False
Now fit the model:
[8]:
# NBVAL_IGNORE_OUTPUT
_ = model.fit(logfreq=200, reg_escape_weight=0.1)
# First fitting site-level model.
# Starting optimization of 1248 parameters at Tue Aug 2 06:37:55 2022.
step time_sec loss fit_loss reg_escape reg_spread reg_activity
0 0.061416 37315 37314 0 0 0.90499
188 11.785 4865.6 4822.7 39.237 0 3.6724
# Successfully finished at Tue Aug 2 06:38:07 2022.
# Starting optimization of 8450 parameters at Tue Aug 2 06:38:07 2022.
step time_sec loss fit_loss reg_escape reg_spread reg_activity
0 0.089303 7299.1 7031.4 264.09 8.7534e-31 3.6724
200 16.369 6821.7 6715.9 80.416 21.021 4.3081
299 24.195 6821.1 6714.3 81.169 21.295 4.3101
# Successfully finished at Tue Aug 2 06:38:31 2022.
Look at output. First we look at the mut_escape_df. The site entries are str:
[9]:
# NBVAL_IGNORE_OUTPUT
assert all(model.mut_escape_df["site"].astype(str) == model.mut_escape_df["site"])
assert set(model.mut_escape_df["site"]).issubset(sites)
model.mut_escape_df.head().round(2)
[9]:
| epitope | site | wildtype | mutant | mutation | escape | times_seen | |
|---|---|---|---|---|---|---|---|
| 0 | 1 | 1 | M | I | M1I | -0.00 | 1 |
| 1 | 1 | 1 | M | T | M1T | -0.01 | 2 |
| 2 | 1 | 2 | F | I | F2I | 0.00 | 2 |
| 3 | 1 | 2 | F | L | F2L | 0.01 | 14 |
| 4 | 1 | 2 | F | S | F2S | -0.01 | 14 |
The same is true for mut_escape_site_summary_df:
[10]:
# NBVAL_IGNORE_OUTPUT
assert all(
model.mut_escape_site_summary_df()["site"].astype(str)
== model.mut_escape_site_summary_df()["site"]
)
assert set(model.mut_escape_site_summary_df()["site"]).issubset(sites)
model.mut_escape_site_summary_df().head().round(2)
[10]:
| epitope | site | wildtype | mean | total positive | max | min | total negative | n mutations | |
|---|---|---|---|---|---|---|---|---|---|
| 0 | 1 | 1 | M | -0.00 | 0.00 | -0.00 | -0.01 | -0.01 | 2 |
| 1 | 1 | 2 | F | 0.01 | 0.03 | 0.02 | -0.01 | -0.01 | 4 |
| 2 | 1 | 3 | V | 0.06 | 0.64 | 0.56 | -0.12 | -0.14 | 8 |
| 3 | 1 | 4 | F | 0.22 | 1.12 | 1.02 | -0.01 | -0.01 | 5 |
| 4 | 1 | 5 | L | 0.02 | 2.08 | 0.94 | -0.89 | -1.78 | 20 |
Now plot the mutation-escape values:
[11]:
# NBVAL_SKIP
# NBVAL_IGNORE_OUTPUT
# we ignore this cell with `nbval` because for some reason it hangs on Travis
model.mut_escape_heatmap(init_min_times_seen=3)
[11]:
And plot the site-escape values:
[12]:
# NBVAL_IGNORE_OUTPUT
model.mut_escape_lineplot(min_times_seen=3)
[12]:
Confirm same results for reference and sequential integer number¶
To demonstrate how the results are the same regardless of which numbering scheme is used for the fitting, we also fit a model with the sequentially numbered variants. This section of the notebook can almost be considered a test rather than an example.
[13]:
# NBVAL_IGNORE_OUTPUT
model_sequential = polyclonal.Polyclonal(
data_to_fit=(
data_to_fit.drop(columns="aa_substitutions").rename(
columns={
"antibody_concentration": "concentration",
"aa_substitutions_sequential": "aa_substitutions",
}
)
),
n_epitopes=1,
alphabet=polyclonal.AAS_WITHSTOP_WITHGAP,
)
assert model_sequential.sequential_integer_sites is True
_ = model_sequential.fit(logfreq=200, reg_escape_weight=0.1)
# First fitting site-level model.
# Starting optimization of 1248 parameters at Tue Aug 2 06:38:39 2022.
step time_sec loss fit_loss reg_escape reg_spread reg_activity
0 0.055258 37315 37314 0 0 0.90499
188 12.714 4865.6 4822.6 39.396 0 3.6721
# Successfully finished at Tue Aug 2 06:38:52 2022.
# Starting optimization of 8450 parameters at Tue Aug 2 06:38:52 2022.
step time_sec loss fit_loss reg_escape reg_spread reg_activity
0 0.074084 7300 7031.3 265.04 8.9894e-31 3.6721
200 17.358 6821.6 6715.8 80.443 21.032 4.3093
245 21.341 6821.4 6715.2 80.692 21.171 4.3108
# Successfully finished at Tue Aug 2 06:39:14 2022.
Now make sure the fitting gives nearly the same result regardless of which site numbering is used. First check the activity values:
[14]:
pd.testing.assert_frame_equal(
model_sequential.activity_wt_df,
model.activity_wt_df,
atol=0.01,
)
Now compare the mutation-escape values. In order to do this, we have to re-number the sequential values to reference numbering:
[15]:
min_times_seen = 50
mut_escape = model.mut_escape_df.drop(columns="mutation")
mut_escape_sequential = model_sequential.mut_escape_df.assign(
site=lambda x: x["site"].map(
site_numbering_map.set_index("sequential_site")["reference_site"].to_dict()
)
).drop(columns="mutation")
# have to use fairly big atol to test this, so also do correlations
pd.testing.assert_frame_equal(
mut_escape,
mut_escape_sequential,
atol=1.5,
)
assert 0.99 < mut_escape["escape"].corr(mut_escape_sequential["escape"])
Test of averaging¶
Just as a toy test, average the reference-site model with itself and make sure it still displays reference-based site numbers:
[16]:
avg_model = polyclonal.PolyclonalAverage(
pd.DataFrame({"number": [1, 2], "model": [model, model]}),
)
[17]:
# NBVAL_IGNORE_OUTPUT
avg_model.mut_escape_heatmap(init_min_times_seen=3)
[17]:
[18]:
# NBVAL_IGNORE_OUTPUT
avg_model.mut_escape_lineplot()
[18]:
You cannot average models that don’t have the same site-numbering setup in terms of sequential or non-sequential:
[19]:
# NBVAL_RAISES_EXCEPTION
# this will give an error
polyclonal.PolyclonalAverage(
pd.DataFrame({"number": [1, 2], "model": [model_sequential, model]}),
)
---------------------------------------------------------------------------
PolyclonalHarmonizeError Traceback (most recent call last)
/tmp/ipykernel_53892/345686643.py in <module>
2
3 # this will give an error
----> 4 polyclonal.PolyclonalAverage(
5 pd.DataFrame({"number": [1, 2], "model": [model_sequential, model]}),
6 )
~/polyclonal/polyclonal/polyclonal_collection.py in __init__(self, models_df, harmonize_to, default_avg_to_plot)
874 harmonize_to = models_df.iloc[0]["model"]
875
--> 876 models_df["model"] = [
877 m.epitope_harmonized_model(harmonize_to)[0] for m in models_df["model"]
878 ]
~/polyclonal/polyclonal/polyclonal_collection.py in <listcomp>(.0)
875
876 models_df["model"] = [
--> 877 m.epitope_harmonized_model(harmonize_to)[0] for m in models_df["model"]
878 ]
879
~/polyclonal/polyclonal/polyclonal.py in epitope_harmonized_model(self, ref_poly)
2160
2161 if self.sequential_integer_sites != ref_poly.sequential_integer_sites:
-> 2162 raise PolyclonalHarmonizeError(
2163 "Models being harmonized don't have same `sequential_integer_sites`"
2164 )
PolyclonalHarmonizeError: Models being harmonized don't have same `sequential_integer_sites`